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ABSTRACT 

An explicit finite difference iteration scheme is developed to 
study harmonic sound propagation in aircraft engine nacelles. To 
reduce storage requirements for large 3D problems, the time depen- 
dent potential form of the acoustic wave equation is used. To insure 
that the finite difference scheme is both explicit and stable, time is 
introduced into the Fourier transformed (steady-state) acoustic 
potential field as a parameter. Under a suitable transformation, the 
time dependent governing equation in frequency space is simplified 
to yield a parabolic partial differential equation, which is then 
marched through time to attain the steady-state solution. The input 
to the system is the amplitude of an incident harmonic sound source 
entering a quiescent duct at the input boundary, with standard 
impedance boundary conditions on the duct walls and duct exit. The 
introduction of the time parameter eliminates the large matrix 
storage requirements normally associated with frequency domain 
solutions, and time marching attains the steady-state quickly enough 
to make the method favorable when compared to frequency domain 
methods. For validation, this transient-frequency domain method is 
applied to sound propagation in a 2D hard wall duct with plug flow. 

INTRODUCTION 

Both steady-state (frequency domain) and transient (time do- 
main) finite difference and finite element techniques have been 
developed to study sound propagation in aircraft nacelles (Baumeister, 
1980a and b and Baumeister and Horowitz 1984). Sound propaga- 
tion with axial variations in duct geometry, mean flow Mach number 
and wall sound absorbers have been considered ( Astley and Eversman, 
1981). To date, the numerical solutions have generally been limited 
to moderate frequency sound and mean flow Mach numbers in two 
dimensional axisymmetric nacelles. Wavelength resolution prob- 
lems have prevented a broader range of applications of the numeri- 
cal methods. A fine grid is required to resolve the short wavelengths 
associated with high frequency sound propagation with high inlet 


Mach numbers. Thus, application of numerical techniques to high 
frequency sound propagation in 3-D engine nacelles has yet to be 
attempted. 

Generally, the number of grid points along the centerline of an 
aircraft nacelle is directly proportional to the sound frequency and 
the nacelle length and inversely proportional to one minus the mean 
flow Mach number. In addition, the number of grid points in the 
transverse direction depends on the radial and spinning mode con- 
tent of the sound source. These dependencies severely limit the 
application of numerical techniques. 

The present paper is the first step in a larger research effort to 
develop efficient numerical techniques to predict high frequency 
sound propagation around 3D aircraft inlet nacelles with large 
subsonic inlet Mach numbers. The paper begins with a description 
of the problem, an explanation of why the transient approach is 
employed, and a description of the grid system and governing 
equations. The bulk of the paper describes the development of a 
stable, explicit finite difference scheme by a transformation of the 
governing hyperbolic wave equation. The scheme is iterated in time 
to converge to the steady-state solution associated with a Fourier 
transform solution. 

NOMENCLATURE 

C dimensionless speed of sound, C^/C 0 ^, eq. (2) 
c steady speed of sound, eq. (1 1) 
c' fluctuating speed of sound, eq. ( 1 2) 

D# dimensional duct height or diameter 
D duct height, D = 1 

d parameter, eq. (58) 
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F source amplitude at duct entrance, eq. (40) 
dimensional frequency 
f dimensionless frequency, f # D # /C c> # 
g parameter, eq. (58) 
h parameter, eq. (58) 

i V-f 

L length of duct, L # /D # , Figure 2 

Mf Mach number at duct entrance 

n unit outward normal 

P* dimensional pressure 

P dimensionless fluid pressure, P # /p 0 # C 0 #2 

p steady fluid pressure 

P' acoustic pressure fluctuation, eq. (19) 

t dimensionless time, f*t # 

tj total dimensionless calculation time 

At time step 

u n ' # normal acoustic velocity, eq. (43) 
x dimensionless axial coordinate, x # /D # 

Ax axial grid spacing 

y dimensionless transverse coordinate, y # /D # 

Ay transverse grid spacing 
Z # impedance 
y ratio of specific heats 
£ dimensionless impedance, Z # /p 0 # C 0 # 
p dimensionless fluid density, p # /p 0 # 
p steady fluid density, eq. (2 1 ) 
p' fluctuating acoustic density 

O dimensionless time dependent flow potential, <J> # /C 0 # D # 


steady mean flow potential, eq. (5) 

<(>' transient acoustic potential, eq. (5) 

<(> transient acoustic potential in frequency space, eq. (30) 

H* spatial potential, eq. (28) 
to dimensionless frequency, 2icf 
V D # V # 

Subscripts 

i axial index, see Figure 2 
j transverse index, see Figure 2 

o ambient or reference condition 

Superscript 

* dimensional quantity 

k time step 

PROBLEM STATEMENT 

The problem under consideration here is the steady-state propa- 
gation of sound, represented by the perturbation acoustic potential, 
through a 2D rectangular duct. The source, noise emanating from fan 
blades in a jet engine inlet nozzle, is represented by specifying the 
pressure distribution at the fan face. The goal of the paper is to 
develop a stable, explicit finite difference scheme that incorporates 
the far field impedance condition applied at the duct exit and rigid 
body boundary conditions on the duct wails. The method is designed 
with the intention of extending the current 2D duct formulation to 
general 3D nacelle design problems with a variety of possible 
boundary conditions in the near and far fields. 

TRANSIENT APPROACH 

In frequency domain approaches, pressure and acoustic velocity 
are assumed to be harmonic functions of time. The matrices associ- 
ated with numerical solutions of the governing equations in the 
frequency domain have extremely large storage requirements. In 
similar 3D electromagnetic applications, Taflove (1991) reports that 
frequency domain approaches are limited to several hundred thou- 
sand unknowns (still a small problem in 3D). Larger problems 
encounter roundoff errors and conditioning problems that prevent a 
reliable solution. 

On the other hand, an explicit transient method generates no 
matrices, and is generally faster than a frequency domain approach 
(Baumeister, 1980a). Miller has developed explicit relationships 
showing the time advantage of transient over steady state techniques 
(Miller, 1988, section V. Some Practical Considerations, C. Discus- 
sion). Consequently, the method of choice in the present paper is a 
time dependent method. 
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GRID SYSTEM 

Three possible mesh systems could be employed to model the 
time dependent acoustic field in an aircraft nacelle; (l)CFD body 
fitted grids, (2) unstructured finite element grids and (3) almost 
highly structured grids. Both the body fitted and finite element grids 
require extensive storage of information about the grid structure. For 
3D problems, this storage requirement places a large burden on 
available storage and speed of operations. Therefore, the almost 
highly structure grid is preferred here. 

The almost highly structured grid is rectangular away from the 
geometry of interest, with nonstandard split or stretched finite- 
volume cells only at the curved surface. Taflove (1991) reports that 
the processing speed is up to 18 times faster than codes using CFD 
body-fitted meshes. Furthermore, artifacts due to refraction and 
reflection of numerical waves propagating across global mesh 
distortions are reduced. 

Figure 1 illustrates an almost structured grid model of an aircraft 
engine nacelle. Only in the vicinity of the center body or the curved 
nacelle tip is the grid nonuniform. For the rectangular duct examples 
considered here, the highly structured grid system shown in 
Figure 2 is employed. 

GOVERNING EQUATIONS 

The governing differential equations for studying acoustic propa- 
gation in inlet nacelles can be formulated in terms of the constitutive 
continuity and momentum equations (Baumeister, 1979) or in terms 
of potential flow (Sigman et. al., 1978). The constitutive equation 
approach can handle a general 3D sheared flow while the potential 
approach is limited to 3D in viscid flow. The major advantage of the 
potential flow approach over the constitutive equation approach 
comes in decreased storage requirements associated with only one 
dependent variable. 

Fortunately, acoustic propagation in inlet nacelles can be reason- 
ably modeled by an inviscid approximation. For single mode JT15D 
engine data, a previous finite element study (Baumeister and 
Horowitz, 1984) employing the potential formulation in the fre- 
quency domain showed good agreement with experimental data - in 
the far field radiation pattern as well as suppressor attenuation. Due 
to this success, the problem under consideration here is formulated 
in terms of an acoustic potential. 

For inviscid, non heat conducting and irrotational flow, the non- 
linear partial differential equation for the flow potential can be 
written in dimensionless form as (Thompson, 1972, pg. 257, 
eq. 5.24) 

f 2 «l> n + f (VO • VO) t + F Vo • V(VO VO) = C 2 V 2 0 (1) 

where 

C 2 =1— (y — l)ffO t +l-VOVoj (2) 

The symbol <5> represents the time dependent potential of the total 
flow field. The speed of propagation of a disturbance is represented 
by C and the frequency of an acoustic source by f. Subscripts 
indicate partial differentiation with respect to subscripted variables. 

The conventional normalization factors used to develop these 


nondimensional equations are given in the NOMENCLATURE. 
However, the normalization of time deserves some special com- 
ment. A common choice for normalizing time is t = C 0 # t # /D # The 
superscript # designates a dimensional quantity while the subscript 
o indicates an arbitrary reference value. With this choice, the 
dimensionless frequency f would not appear in eqs. (1) or (2). 
However, in this paper, the dimensional frequency f* of the forcing 
acoustic signal was chosen to normalize time, so that t = f*t # . As a 
result, the time t indicates the number of complete acoustic cycles 
that have occurred since the start of the solution process. This is 
advantageous because the total time of the numerical calculation can 
generally be set independently of the frequency of the acoustic 
signal. 

Rewriting eqs. (1) and (2) in two dimensional rectangular coor- 
dinates yields: 

f 2 <D„ + f(20 x ® xt + 20 y O yt ) + <D 2 <D XX + 2<D x O y <I> xy 

2 21 \ (3) 
+ <pJ<l> yy =C 2 (<I> xx +O yy ) 

C 2 =l-(y-l)(«*t+^+^) (4) 

To obtain the acoustic solution, the flow potential is 
rewritten as the sum of a steady mean flow potential O (x,y) and an 
acoustic potential <|)'(x,y,t); that is 

4>(x,y,t) = <D(x,y) + <t>'(x,y,t) (5) 

To simplify the algebra in the linearization of eq. (3), first eq. (4) is 
linearized and the steady mean flow equation is developed. Substi- 
tuting eq. (5) into eq.(4) and neglecting second order terms, the 
speed of the disturbance can be expressed in terms of the steady and 
fluctuating potentials as 



which can be written in compact form as 


C 2 = B 2 +B' 2 

(7) 

where 


B 2 =l-i(Y-l)(^x+^y) 

(8) 

B '2 = _( Y _iX f ^;+o x 

(9) 

To linearize the speed of sound, define 



C = c + c' = ^B 2 +B' 2 =B 1 + E = b + ~ (10) 

\ B z 2 B 
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Therefore, the mean speed of sound depends on the velocity field 
and can be expressed as 

c = B = { 1 -i-( T - 1 )(oJ +oj)} (11) 

while the perturbation of the sonic velocity depends both on the 
mean flow field and the acoustic field as 


Assuming that 

P(x, y, t) - P(x, y) + P'(x, y, t), (19) 

p(x,y,t)=p(x,y) + p'(x,y,t), (20) 

substituting eqs.(5), (19), and (20) into eq. (18) and neglecting 
higher order fluctuation terms yields for the acoustic pressure P / 


c ' = -(¥)( f ^ + ^ + Vy) < 12 > 

The differential equation describing the mean flow can now be 
expressed in terms of the mean sonic velocity. Substituting <l> into 
both eqs. (3) and (4) yields 

+2^ y O xy + 0 2 0 yy - c 2 (o xx +® yy ) = 0 (13) 

Now, the linearized wave'equation can be conveniently estab- 
lished. Substituting eqs. (5) and (7) into eq.(3), factoring out the 
steady contribution from eq. (13), and neglecting nonlinear products 
of the acoustic quantities yields 

0 = f : Vtt - (c 2 - «?)♦'„ - (c 2 - O 2 ) 4>' yy + 2 O x 0 y *' xy 

+ 2f O x fxt + 2f O y 4.^ t + 2(o x O xx + O y O xy )<t> x 

+ 2(0 X 0 X y + Oy^yy )^ - (Y — 1) 

X (ffj+O.^+Oyljj'y) (Oxx+^yy) (14) 

For the special case of plug flow, the mean flow terms in eq. (14) 
become 

*y=*xy=*xx=0 *x=M f (15) 


= P,' = -fft<-(‘M>;+<Vy) t ( 21 > 

Finally, at the fan face (x = 0), assume 

O y =0 at x = 0 (22) 

O x =M f at x = 0 (23) 

and integrate with respect to time to obtain 

V(0,y,t) = -f<t>;-M f <|>; (24) 

p 

Eq. (14) and eq. (24) are the basic equations used to establish a 
finite difference formulation for the rectangular duct shown in 
Figure 2. However, care must be taken when discretizing derivative 
terms to insure that the resulting scheme is both stable and explicit. 
Note that it is easy to develop a stable implicit method, but this would 
yield a matrix formulation that is no better than a frequency domain 
approach. It turns out that the proper treatment of the mixed time and 
space derivative terms (which appear on the second line of eq. (14)) 
is critical for maintaining stability. 

The implicit formulation is obtained by approximating the mixed 
partials by (Baumeister, 1980b, eq. (16), and Abramowitz and 
Stegun, 1964, pg 884, eq. 25.3.27) 


Substituting eq. (15) into eq. (14) yields 


t>» k j +1 + » i-l, j + ♦& . j + 1 - 20 i ft - 


0 = f 2 - (c2 - M 2 y xx - c2 *' yy + 2 f M f 4>' xt (16) 


+ ♦vj 1 + ’t’llj+l + 4>'iij-l - 2< I>|J - - •t’jjy 


c 2 =1--(y-1) M 2 


Generally, the pressure at the fan face is specified. Therefore, the 
relationship between pressure and the potential needs to be deter- 
mined from the conservation of momentum. The dimensionless 
form of the conservation of momentum equation can be written as 


where i and j denote the space indices for the nodal system shown in 
Figure 2, k is the time index defined by 


t k+1 =t k +At 



VO), 


(18) 


and Ax, Ay, and At are the space and time mesh spacings 
respectively. 
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Baumeister( 1980b) showed that eq. (25) can be used in an explicit 
fashion for ID plug flow problems in a duct. However, if the flow 
is not one dimensional or if the region exterior to the duct is included, 
the scheme must be implicit. This approach is inappropriate for 
general 3D problems. Fortunately, this problem can be circum- 
vented by transforming the potential and consequently modifying 
the governing equation. The details follow in the next section. 

TRANSFORMATION TO FREQUENCY SPACE 

There are several ways to develop a frequency domain formula- 
tion for the general 2D acoustic wave eq. (14) or the plug flow 
simplification eq. (16). The Fourier Transform can be applied if the 
potential has a multi-frequency content. In the monochromatic case 
(Temkin, pg 52, eq. 2.5.1), this is equivalent to assuming that 

<t>'(x,y,t) = 4'(x,y)e- i0)V = 4'(x,y)e~ i2m (28) 

which, in the case of plug flow (from eq. (16)), yields 

0 = |c 2 -Mf j'F xx +c 2K P yy +(o 2 'J / + i2coM f 4 / x (29) 

This equation would be solved numerically using a linear system 
of equations. However, the associated matrix is not positive definite, 
which can lead to numerical difficulties, and which preclude the use 
of iterative techniques (see TRANSIENT APPROACH). Therefore, 
it is desirable to develop an explicit finite difference scheme to avoid 
the use of matrices. The situation is complicated by the fact that 
boundary conditions can introduce instabilities (Baumeister, 1982 
and Cabelli, 1982) in the solution process. In time-dependent form, 
eq. (14) or (16) cannot easily be discretized in such a way that the 
resulting finite difference scheme is both stable and explicit in the 
presence of flow (it is possible to obtain reasonable results in the no- 
flow case). 

The resolution of these difficulties is achieved by modifying the 
monochromatic transformation (28) to 

<|>'(x,y,t) = <t>(x,y,t)e -lt0 1 = <|>(x,y,t)e -l2,n (30) 

This differs from the classical monochromatic transformation in that 
the amplitude <|) (no prime) is no longer assumed to be independent 
of time, as in eq. (28). Employing eq. (30), the time derivatives in 
eqs. (14) and (16) can replaced by the following relationships: 

<)>{ =[-i27i<t> + <t> t ]e _l27Ct (31) 

f xt =[-i2^ x+ 4> xt ]e- i2ro (32) 

= -2i2m() t -(2rt) 2 4»)e 12 ” (33) 

Under this transformation, the general equation (14) and the plug 
flow equation ( 1 6) become, respectively 


f 2 <1>„ -[i2fo) + (y-l)f(<I>„ +«> yy )]<t>, +2ra>A, 

+ 2fO y <|> yl = (c 2 - O 2 )<|> xx +(c2-0 2 )<t> yy 
-2<i> x C> y <t> xy +C0 2 <t> + i20)O x ((> x +i2(DO y <t> y 

+O y O Jty )<(» x -2 (1,0,,+®,^) 

X 4> y +(Y-l)(- iaX t , + ^x't>x +4Vt> y )(«\x + *\y) 

(34) 

f 2 <t> tt -2ifa*t +2fM f <D xt =(c2 -M 2 )* xx 

+ c 2 <(>yy + co 2 <j) + i2coM f <(> x (35) 

To see the relationship between the two transforms (eq. (28) 
and (30)), consider, for simplicity, the case of plug flow. The 
monochromatic transform yields eq. (29), while the transient trans- 
formation yields eq. (35). The only difference is in the presence of 
the time derivative terms on the left hand side of eq. (35). Physically, 
the time dependence in 4>(x,y,t) is caused by assuming that the duct 
is quiescent at time 0, and that the source is turned on at that instant. 
A series of numerical calculations was performed to investigate the 
relationship between y and <{► as time progresses. It was verified that 
<j > converges to the steady state y, so that 

lim <|)(x,y,t) = 'F^y) (36) 

t—*oo 

At present, a formal proof of convergence is not available. 

At present, transient solutions in frequency space are of no 
interest. Therefore, approximations to eqs. (34) or (35) can be made, 
by modifying the time derivative terms, as long as the following two 
factors are taken into account: 1 ) the resulting scheme must be both 
stable and explicit, and 2) there must be at least one time derivative 
term in the equation to ensure that the scheme can be marched 
through time to obtain the steady state solution. The mixed space- 
time derivative term in eq. (35) prevents an explicit finite difference 
representation of the governing equation, and is therefore dropped. 
Furthermore, the hyperbolic time term in eq. (35) is also dropped to 
further simplify the governing equation. For plug flow, this produces 
a parabolic “wave” equation of the form: 

-2ifOK)> t =(c 2 -M 2 j<|> xx + c 2 <|> yy +co 2 <j) + i2coM f <|> x (37) 

This type of differential manipulation is often associated with 
preconditioning of both time dependent partial differential equa- 
tions (Turkel, Fiterman and Leer, 1 994) and time independent partial 
differential equations (Turkel and Amone, 1993) to accelerate 
convergence to the steady state solution. Here, though, the goal is 
obtaining a stable, explicit scheme. 

In the absence of mean flow, the parabolic wave equation takes on 
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the form 

-2ifax|>, =<(>„+ <t>yy + (0 2 <t> (38) 

The right hand side of eq. (38) becomes the Helmholtz equation. In 
the Fourier transform approach, Bayliss, Goldstein, and Turkel 
(1982) have developed an iterative approach to its associated 
matrix; however, the Gaussian elimination approach is still the 
method of choice. In effect, the time iteration procedure presented 
in this paper offers a very simple iterative solution to this classic 
equation. 


INITIAL AND BOUNDARY CONDITIONS 

The duct is assumed to be quiescent at time 0, so that the initial 
condition is 

4>(x,y,0) = 0 (39) 

As the equation is iterated in time, the solution builds up to the steady 
state solution. 

At the duct entrance, (x= 0), the potential is given by 


component of the particle velocity normal to its surface, pointing 
into the medium characterized by Z #> . The velocity u^* is positive 
if its direction points into the second medium, ie., to the outside of 
the surface that contains the incident wave. 

In dimensionless form, eq. (43) can be written as 





p' 

d<(>' 

dx 


(44) 


Substituting in the relationship between pressure and potential for 
plug flow as given in eq. (24) yields 


or 


’ r A. ' 


- f »; 

— + M f 
P 


(45) 


(46) 


<t.'(0,y,t) = F(y)e- i2m (40) 

and through eq. (30) to 

<t>(0,y,t) = F(y) (41) 

If the pressure at x=0 is specified as the boundary condition, then the 
potential is related to the pressure directly through eq. (24). 

The wall and exit boundary conditions in hyperbolic space 
generally require special considerations. Rigorous treatment of time 
dependent boundary conditions for hyperbolic systems is quite 
involved (Thompson, 1987). In time dependent duct propagation 
problems, impedance concepts have generally not proven useful 
(Banks, Propst and Silcox, 1991). However, in the frequency 
formulation, Baumeister & Horowitz (1984) found impedance 
concepts very useful in describing the operation of real jet engine 
noise suppressors. Thus, a major advantage of the transient- fre- 
quency formulation presented here is the capability of using the 
impedance and gradient conditions as developed for the frequency 
domain. 

The hard wall condition is 


Since the equations are being developed for harmonic flow, the 
following relationship can be used to simplify eq. (46): 

<K| t ^„=-i2rotK : - i2m (47) 

Therefore, the relationship between the gradient of $ and the 
impedance in eq. (46) can be written as 

ico 

<t> x =7 * (48) 

± + M f 
P 

For plane wave propagation, the plug flow solution for the 
acoustic potential is 

KD 

<(>'(x,t) = e‘ +M,> e~ i2m (49) 

Substituting eq. (49) into eq. (45) yields 


V<f> n =0 


(42) 


C=pc 


(50) 


where n is the unit outward normal. 

To simulate a non-reflective boundary at the duct exit, the 
difference equation is expressed in terms of an exit impedance. For 
the examples presented in this paper, the duct exit impedance is 
taken to be that for a plane acoustic wave. The impedance Z # is 
defined as (Skudrzyk, eq, (34) section 15.4, pg. 299) 


P t fl 

Z" =— = 


p/tl 




(43 


where P' # is the acoustic pressure at the interface and u n ' # is the 


Substituting eq. (50) into eq. (48) gives the value of the exit gradient 
at the end of the finite difference domain: 


ICO 


c + M f 


(51) 


This gradient is used to simulate the non-reflective exit condition. 

This approach can be used only to simulate plane waves in ducts 
or at a termination deep in the far field from an engine nacelle. For 
multimode propagation in ducts or grid termination in the near field, 
the modal -element method has been found to be very useful. That 
method can be used within this transient- frequency framework. In 
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acoustics, Astley and Eversman (1981) and Baumeister and Kreider 
(1993) applied the modal element method to duct acoustic and 
scattering problems. In conventional CFD, Baumeister and 
Baumeister (1994) also applied the method to potential flow over a 
cylinder within in a duct. 

FINITE DIFFERENCE EQUATIONS 

The finite difference approximations determine the potential at 
the spatial grid points at discrete time steps & = kAt. Starting from the 
known initial conditions at t = 0 and the boundary conditions, the 
algorithm marches the solution out to later times. 

Away from the duct boundaries, as shown by the cell in Figure 2, 
each partial derivative in eq. (37) can be expressed as follows: 


(1980a & 1980b) gives precise details for generating the time 
difference equations at the boundaries. 

STABILITY 

A von Neumann stability analysis (Lapidus and Pinder, 1982) 
indicates that the method is conditionally stable, subject to the 
conservative condition 


2 

7_df +| 

r <n 

2" 


cof 

IaxJ + 

Uy j 


f Ax 


.k+l A k-\ 

♦o -*0 

2At 

(52) 

Ax 2 

(53) 

Ay 2 

(54) 

2Ax 

(55) 

M .1? 

II 

(56) 


Substituting these expressions into eq. (37) yields 



where 

d 2 =c*-M 2 

h = icof (58) 

g = 2icoM f 

Eq. (57) is an explicit two step scheme. At t = 0, field values at 
t* -1 are assumed zero because the initial field is quiescent. 

The expressions for the difference equations at the boundaries are 
complicated somewhat by the impedance conditions. However, a 
simple integration procedure resolves this problem. Baumeister 


In a typical application, CO, f, and Mf are set by the engine operating 
conditions. Next, the grid spacing parameters Ax and Ay are set to 
accurately resolve the estimated spatial harmonic variation of the 
acoustic field. Finally, At is chosen to satisfy eq. (59). 

In the von Neumann analysis, conditional stability means that the 
amplification factor, which describes how errors propagate from one 
time step to the next, has magnitude one. Thus, when inequality (59) 
is satisfied, errors are not magnified or diminished in magnitude. 
This is a desirable property, since the numerical formulation can not 
distinguish between an error and a small acoustic mode. 

The von Neumann stability analysis does not take into account 
boundary conditions. For stability, gradient boundary conditions 
require the use of smaller At than predicted by inequality (eq. (59)). 

NUMERICAL EXAMPLES 

In the three examples that follow, the parabolic transient- frequency 
domain results are compared to the exact results of the steady Fourier 
transformed solutions, given by eq. (49). 

The following problem is considered: a plane wave propagates 
from the left into a quiescent duct of length one, and the acoustic 
potential field is to be computed in the duct. Note that, boundary 
conditions can introduce instabilities (Baumeister, 1 982 and Cabelli, 
1982) into otherwise stable finite difference schemes. Therefore, it 
is important to test the proposed method for convergence in time to 
the steady state solution in the absence of the exit boundary condition 
(eq. (51)), and to test independently the effect of the exit boundary 
condition itself on the solution. 


Semi-Infinite Duct 

In this example, the computational boundary is set far enough 
away from the true boundary x=l that any artifacts arising from 
imperfections in the exit boundary condition do not affect the 
solution in [0,1]. The numerical solution propagates one node per 
time step, so setting the boundary at x=50 with step Ax =0.05 
provides a sufficient number of time steps to gauge the convergence 
of the method before any artifacts might reflect back from the 
computational boundary. 

The numerical and exact results are compared in Figure 3, for no 
flow (Fig. 3a) and for Mach number Mf =-0.5 (Fig. 3b). The 
frequency is normalized to f=l. Both cases show excellent agree- 
ment. The total calculation time was t T =5.0. 
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Finite Duct L=1 

In this example, the computational boundary is moved up to the 
true boundary x=l to examine the effect of the exit boundary 
condition (eq. (5 1 )). The frequency is normalized to 1 . Again, two 
cases were considered - no flow (Mf = 0) and Mach number 
Mf =-0.5. The results are shown in Figure 4a and 4b, respectively. 
The numerical results again match well with the exact results, but it 
is clear that the exit boundary condition does have a slight degrading 
effect on the solution. Notice also that the time step has been 
decreased here, which tends to increase the execution time; how- 
ever, the computational domain is smaller, which tends to decrease 
the execution time. The total calculation time in this example was 
tj = 4.0. 

CONVERGENCE RATE 

In this example, the convergence rate is studied for the region 
x = 0 to x = 1 using the semi-infinite duct. The results are shown in 
Figure 5 for the real and imaginary components of the potential and 
Figure 6 for the magnitude of the potential. As seen in these figures, 
the numerical solution quickly and accurately converges to the exact 
steady state solution. 

For this plane wave propagation, the exact solution to the hyper- 
bolic governing equations predicts the arrival of the “steady state” 
solution when t=l (Baumeister, 1983, eq. (28) renormalized). As 
seen in Figures 5 and 6, the parabolic solutions converges to the 
Fourier transformed results after a time of t=2 has elapsed. Thus, the 
parabolic marching scheme requires more time steps to converge 
than a hyperbolic method. 

For parabolic partial differential equations with real coefficients, 
disturbances propagate with infinite speed (Morse and Feshbach, 
1 95 3, pg. 865). In the present calculations, numerical solutions to the 
parabolic eq. (37) have this trait. The finite difference values of the 
potential propagate throughout the domain at the numerical velocity 
Ax/At rather than the speed of sound. This characteristic may 
account for the slower convergence time for the parabolic formula- 
tion as compared to the hyperbolic formulation of the appropriate 
governing equations. In the parabolic scheme, acoustic energy is 
numerically transferred ahead of the true wave front. Nevertheless, 
the parabolic scheme provides the correct steady state solution using 
minimal computer storage with run times comparable to other 
methods. 

HISTORICAL PERSPECTIVE 

With the transient-frequency approach developed in this paper, 
three different finite difference/element solution techniques are now 
available to solve the hyperbolic wave equation that describes 
acoustic propagation in jet engine nacelles. This is outlined in Figure 
7 for the zero mean flow case. 

The Fourier transform of the wave equation was the first numeri- 
cal approach used to study sound propagation in jet engine ducts 
(Baumeister and Bittner, 1973 and Baumeister and Rice, 1973). This 
steady state approach is outlined on the right side of Figure 7. The 
governing hyperbolic wave equation is transformed to the elliptic 
Helmholtz equation. Finite difference (FD) and finite element (FE) 
numerical formulations have been employed to solve this equation. 
After applications of the boundary conditions (Fig. 7; [BC]), the 
associated finite difference or finite element global matrix is solved 


for the velocity potential (or pressure). Because the matrix form of 
the Helmholtz equation is not positive definite, matrix elimination 
solutions are generally employed. This requires extensive storage, 
as discussed in TRANSIENT APPROACH. Conveniently, the 
steady state approach allows the direct calculation of the potential 
(or pressure) fields. 

In the inlet to a turbojet engine, the dimensionless frequencies f 
can be on the order of 30 to 50 for the higher harmonics of the blade 
passing frequency. The storage requirements and associated com- 
puter run times for these high frequencies make computations 
expensive or even impossible. To make the numerical solutions 
more cost effective, grid saving approximations to the governing 
Helmholtz equation have been used (Fig. 7; [Approximate]). 
Baumeister ( 1 974) employed the wave envelope theory while Hardin 
and Tappert (1973) developed a similar approach for underwater 
sound propagation with the addition of a parabolic (space) approxi- 
mation. Candel (1986) presents an extensive discussion of the 
contemporary research in this area and a detailed development of the 
parabolic equation method (PEM). 

The transient solution to the wave equation was the second 
numerical approach used to study propagation in jet engine ducts, 
which is shown in the left-hand column of Figure 7. To eliminate the 
matrix storage requirements, Baumeister developed time dependent 
finite difference numerical solutions for noise propagation in a two- 
dimensional duct without flow (1980a), with parallel shear flow 
(1979) and with axi symmetric flow (1980b). 

Sound is introduced as a boundary condition at the duct entrance. 
The initial conditions generally assume a quiescent duct. Finite 
difference (FD) approximations to the hyperbolic wave equation are 
then solved by an iteration process. The calculations are run until the 
initial transience dies out and steady harmonic oscillations are 
established. Finally, the transient variable <t>' is transformed into the 
steady state variable y associated with the solution of the Helmholtz 
equation. As with the steady state approach, simplification to the 
governing equations (Fig. 7; [Approximations]) can reduce com- 
puter storage and run times (Baumeister, 1986). 

The third option, the transient-frequency technique, is illustrated 
in the central column of Figure 7. This fully explicit iteration method 
eliminates the large matrix storage requirements of steady state 
techniques and allows the use of conventional impedance condi- 
tions. As time increases, the iteration process directly computes the 
steady state variable \| f. 

CONCLUDING REMARKS 

A transient-frequency domain numerical solution of the potential 
acoustic equation has been developed. The potential form of the 
governing equations has been employed to reduce the number of 
dependent variables and their associated storage requirements. This 
fully explicit iteration method represents a significant advance over 
previous steady state and transient techniques. Time is introduced 
into the steady state formulation to form a hyperbolic equation. A 
parabolic approximation (in time) to this hyperbolic equation is 
employed. The field is iterated in time from an initial value of 0 to 
attain the steady state solution. 

The method eliminates the large matrix storage requirements of 
steady state techniques in the frequency domain but still allows the 
use of conventional impedance conditions. Most importantly, the 
formulation is fully explicit under flow conditions. In each example 
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provided, the numerical solution quickly and accurately converges 
to the exact steady state solution. 
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Far field boundary 



Far field boundary 


Figure 1. — Almost-completely structured FD-TD mesh for air- 
craft inlet acoustic nacelle model. 



Figure 2. — Structured FD-TD mesh for rectangular duct. 


Real-imag, <J> and i|# Real-imag, 4> and i|i 
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X axial coordinate 
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X axial coordinate 

Figure 3. — Analytical and numerical potential profiles along 
wall for plane wave propagating in a semi-infinite hard 
wall duct (f = 1). (a) M f = 0 (Ax = 0.05, At = 0.003, t T = 5.0). 
(b) M f = -0.5 (Ax = 0.025, At = 0.001 , t T = 5.0). 
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Real-imag, 4> and Real-imag, <J> and i \t 




Figure 4. — Analytical and numerical potential profiles along 
wail for plane wave propagating in a hard wall duct of unit 
length and a non -reflecting exit condition (f = 1). (a) M f = 0 
(Ax = 0.05, At = 0.0001 , t T = 4.0). (b) M f = -0.5 (Ax - 0.025, 
At = 0.00001, t T = 4.0). 
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Figure 5. — Developing history of disturbance propagation 
in Fourier transformed domain as a function of time. 

(M f = 0, Ax = 0.05, At = 0.003). (a) t = 0.25. (b) t = 0.5. 

(c) t = 0.75. (d) t = 1.00. (e) t = 1 .5. (f) t = 2.0. (g) t = 2.5. 
(h) t = 3.0. 
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Figure 6. — Developing history of magnitude of disturbance 
propagation in Fourier transformed domain as a function of 
time. (Mf = 0, Ax = 0.05, At = 0.003). (a) t = 0.25. (b) t = 0.5. 
(c) t = 0.75. (d) t = 1 .00. (e) t = 1 .5. (f) t = 2.0. (g) t = 2.5. 

(h) t = 3.0. 
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Figure 1 — Alternate finite difference/element methods in solving wave equations. 
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